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ORIGINAL ARTICLE 

Predicting cooperative drug effects through the quantitative 
cellular profiling of response to individual drugs 

J Zhao 1 , X-S Zhang 1 and S Zhang 1 

Quantitative prediction of cellular responses to drugs and drug combinations is a challenging and valuable topic in pharmaceutical 
research. In the past decade, microarray technology has become a routine tool for monitoring genome-wide expression changes 
and has been widely adopted for exploring drug response in the pharmaceutical field. However, how to predict the synergistic 
effect of drug combinations using microarray data is a challenging task. In this article, we report a simple prediction framework 
based on the genome-wide and quantitative profiling of cellular responses to individual drugs. By exploring the differential 
expression profiles, our correlation-based strategy can reveal the synergistic effects of drug combinations. The comparison 
with gold-standard experimental results demonstrates the strengths and weaknesses in relation to prediction based only on 
cellular response to individual drugs. Specifically, the prediction strategy may work for a drug combination whose individual 
drugs show related transcriptomic mechanisms but not for others. 
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An understanding of the complex biological responses of 
the human organism to drugs is key to investigating the effi- 
cacy and safety of compounds in drug development. Many 
genomic features, including DNA methylation patterns, mes- 
senger RNA levels, and protein expression or metabolite 
profiles, may be used for monitoring biological responses. 
Microarray is currently the least expensive high-throughput 
technology for simultaneously monitoring genome-wide 
expression profiling of transcriptional response to drugs. 

Microarray data have been systematically explored in 
model organisms to elucidate the drug mechanism of action, 
and coexpression analysis enables the inference of func- 
tional roles for genes that respond coherently to drug pertur- 
bations. The National Cancer Institute's NCI-60 project and 
the Connectivity Map have extended the concept of genome- 
wide gene expression profiles of drug response to human 
cell lines. 12 The NCI-60 project screened 60 human tumor 
cell lines against more than 100,000 compounds and con- 
structed a public repository for basal gene expression and 
drug sensitivity information. The Connectivity Map project 
generated genome-wide expression profiles both before and 
after drug treatment for 1 ,309 compounds and constructed 
a drug network by comparing ranked lists of up- and down- 
regulated genes. 3 More recently, Iskar ef a/. 4 identified a large 
set of drug-induced transcriptional modules from genome- 
wide microarray data of drug-treated human cell lines and 
rat liver and examined their conservation across tissue types 
and organisms. The large-scale Cancer Cell Line Encyclo- 
pedia project extended such analysis and generated gene 
expression, chromosomal copy number, and massively 
parallel sequencing data from 947 human cancer cell lines. 
They analyzed the data using pharmacological profiles for 
24 anticancer drugs across 479 of the cell lines and identi- 
fied genetic-, lineage-, and gene expression-based predic- 
tors of drug sensitivity. 5 All these efforts have enabled the 



systematic translation of molecular data into biological knowl- 
edge and therapeutic possibilities. 

However, it is well known that redundancy and multifunc- 
tionality are inherent characteristics of biological systems. 
Most complex diseases are usually caused by the inter- 
weaving of multiple biological processes that are redundant 
and robust to the perturbation of one single drug. 67 Thus, 
selecting proper drug combinations that can yield a syner- 
gistic effect is thought to be an effective way of countering 
biological buffering and increasing therapeutically relevant 
selectivity. 8-12 By modulating the activity of distinct proteins 
through drug combinations, we can gain not only the poten- 
tial of overcoming the redundancy underlying pathogenic 
processes but also fewer side effects. 1314 For example, estro- 
gen, which modulates feeding behavior and energy expendi- 
ture, has been suggested as a therapy for obesity and type 2 
diabetes; however, it can also cause unwanted effects, such 
as tumor growth promotion. By chemically linking estrogen 
to the hormone GLP-1 to selectively target metabolically 
relevant tissue, Finan ef a/. 15 obtained a more selective and 
potent metabolic drug with fewer side effects. 

Despite these recent advances, understanding the molec- 
ular mechanisms underlying drug combinations is challeng- 
ing. High-throughput screening is useful to identify possible 
drug combinations; however, it is impractical to screen all 
possible drug combinations for all possible indications. We 
still have no practical experiment strategy to identify possible 
effective drug combinations. Several computational methods 
have recently been proposed. 16-22 Molecular and pharmaco- 
logical features of drugs are the common basis of these com- 
putational methods. By integrating several of these features, 
such as their targets, indications, or chemical structures, 
researchers can identify the similarity between each drug 
pair. The drug pairs with significant high similarity can be con- 
sidered candidates for possible effective drug combinations. 
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Signaling pathways and human protein-protein interaction 
networks are two other resources on which many models 
have been based. Jin et a/. 18 proposed an enhanced Petri- 
Net model to identify the targeted signaling networks of two 
drugs separately and their combination. By comparing identi- 
fied drug effects, they could verify whether a drug combina- 
tion shows synergistic effect. 

Only a few studies have been devoted to predicting syn- 
ergistic effects from gene microarray data. 1819 Their com- 
mon strategy is comparing the drug effects derived from the 
microarray treated by drug combination with those from the 
corresponding two microarrays treated by two drugs sepa- 
rately. In our study, we investigated whether cellular responses 
of drug combinations can be modeled quantitatively by the 
cellular profiling of response to individual drugs. We demon- 
strated quantitative modeling based on the expression profiles 
acquired after treatment with 1 4 individual compounds from the 
NCI-DREAM Drug Sensitivity Prediction Challenge. Because 
the challenge is greater in the absence of cellular response 
profiling of drug combinations, we needed to find common 



features from their respective response data. Our hypothesis 
underlying the current method is that the differential expressed 
genes after administration of one drug can express the action 
mechanism of this compound. The combination of two drugs 
with correlated or anticorrelated differential expression profiles 
may have, respectively, a synergistic or an antagonistic effect. 

RESULT 

Overview of the analysis 

Recently, the Dialogue on Reverse Engineering Assessment 
and Methods (DREAM) (http://www.the-dream-project.org/ 
challenges/nci-dream-drug-sensitivity-prediction-challenge) 
consortium designed a community-based, collaborative com- 
petition for the systematic and objective validation of com- 
putational methods to predict the cooperative effects of 14 
distinct drugs/compounds (M. Bansal etal., personal commu- 
nication). This project generated gene expression data only 
for samples treated by individual drugs, not pairwise combi- 
nation in a human j3-cell lymphoma cell line lymphoma cell 
line (Figure 1a), which is a key difficulty. Predictions from 31 
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Figure 1 (a) Illustration of the data and the research goal of this study. The study was based on experiments that yielded molecular information 
and sensitivity data for a diffuse large B-cell lymphoma cell line (LY3) after the application of 14 different therapeutic drugs, profiled at three 
different time points after drug administration, and baseline expression in growth media. All gene expression profiles of drug-treated cells 
are provided in triplicate. We used the provided gene expression profiles to predict the order of efficacy of all pairs of drug combinations from 
the most synergistic to the most antagonistic, (b) The flowchart of our method used for predicting and ranking cooperative effects of drug 
combinations. After preprocessing, we extracted signature genes to reflect the effects of individual drugs and then calculated the expression 
correlations between their overlapping signature genes to capture the effects of drug combinations. 
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different methods have been experimentally assessed, and 
the present method (with slight modifications here) has been 
ranked fourth. Such results confirm that our method is a simple 
but promising one to predict cooperative effects of drug pairs. 

The goal of this study was to rank the activity of all 91 
pairwise combinations of 14 distinct drugs/compounds in 
a human diffuse large B-cell lymphoma (DLBCL) cell line 
(LY3) from the most synergistic to the most antagonistic 
(Figure 1a). We used only the baseline expression profile 
and the expression profiles treated with each individual drug/ 
compound, with a limited number of conditions and replicates 
to predict synergistic effects of drug/compound combinations. 
The basic assumption of our method is that the expression 
properties of differential expressed genes obtained by com- 
paring the expression profile of a cell line before and after 
treatment with a drug reflect the effect of this drug to that 
cell line. We further assume that the correlation of expres- 
sion profiles of differential expressed gene sets for two drugs 
will reflect the synergistic and antagonistic effects of the two 
drugs (see Figure 1b and Methods). All results and codes 
are provided in the Supplementary Information. 

The predictions were then evaluated against an experimen- 
tally assessed "gold standard" generated by evaluating com- 
pound combination activity in vitro (M. Bansal et al., personal 
communication). A measure of the concordance between the 
predicted rank list and the gold standard known as the con- 
cordance index (c-index) was computed to score the predic- 
tion. In essence, the c-index computes the fraction of item 
pairs whose relative order in the two rank lists is coincident. 

Differential expressed genes indeed provide valuable 
information for predicting cooperative drug effects 

In our method, we set a P value threshold of 0.05 to choose 
the significantly differentially expressed genes as signature 
genes used for the following correlation analysis. Using this 
signature gene set, we can get a c-index of 0.599, which 
is very competitive compared to the top two methods, with 
c-indexes of 0.613 and 0.605, respectively (M. Bansal era/., 
personal communication). To show its effectiveness, we also 
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Figure 2 (a) The statistical significance of random perturbation 
experiment of the chosen differential expressed genes. The 
histogram is plotted based on the concordance index (c-index) 
values of 1 ,000 randomly selected signature genes. The red line 
represents our results based on the chosen differential expressed 
genes. This result is significantly larger than that of randomly 
selected gene signatures with P value = 3.39 x10~ 5 (Mest). (b) The 
c-index values of different differential expressed gene sets based on 
the P value thresholds, confirming that P value = 0.05 is one of the 
"optimal" selections. 



adopted a pathway-based approach for the current task. Spe- 
cifically, we mapped the differential score of expressed genes 
using the Human Protein Reference Database (HPRD) and 
then detected the significant "active" pathway which corre- 
sponds to the maximal weighted connected subnetwork by 
the BioNet method. However, it only shows equivalent per- 
formance compared to the correlation-based strategy for 
this task the in the current setting with weighted c-indexes of 
0.608 and 0.599, respectively. 

We conducted a random perturbation experiment to test 
the significance of the selected differentially expressed 
genes. For each drug combination, we randomly extracted the 
same number of signature genes as the overlap we obtained 
using a P value of 0.05. Then we calculated the correlation 
of expression profiles of these two drugs based on the ran- 
domly selected signature genes. The selection and computa- 
tion were repeated three times, and the one with the maximal 
absolute value was selected as the synergistic score. Then 
we used the scoring scripts to calculate one c-index. The 
entire procedure was repeated 1 ,000 times (the distribution 
of these c-index values is plotted in Figure 2a). We can see 
that the differential analysis indeed has the power to iden- 
tify the most relevant signature genes for the predicting task 
(P value = 3.39 x 10 5 ; ?-test). 

We further assessed the robustness of the P-value thresh- 
old to see how it affected the results. Figure 2b shows the 
values of the c-index for each threshold. We obtained the 
more informative predictions with a P value = 0.05 than with 
others, especially those with a larger value. We should also 
note that with a P value <0.001, we could get only a limited 
number of differential genes that could not be used for the 
correlation analysis. 



Table 1 The top-ranked drug combinations per the gold-standard method and 
their predicted ranking by our method 
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"Shows the modified version. 
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The drug combinations with top consistent predictions 
show similar functional mechanisms 

We further inspected our predictions compared with the 
experimental "gold-standard" rank (Table 1). We found that 
some top predictions that also have relative top ranking by 
the gold standard show highly similar functional mechanisms. 
For example, the predicted top 1 drug combination (mitomy- 
cin and etoposide) was ranked 3 in the gold-standard list. 
We analyzed their overlapping signature genes by functional 
enrichment analysis and extracted connected subnetworks 



by mapping them onto a protein interaction network. We found 
that the common signature genes are enriched in DNA repair, 
DNA metabolic process, DNA replication, cell cycle biological 
processes, and the p53 signaling pathway (Figure 3). More 
specifically, we observed that treatment with these two drugs 
can cause upregulation of p53 signaling genes (CDKN1A, 
PPM1D, TP53I3, SERPINE1, DDB2, MDM2, RRM2B, FAS, 
GADD45B, SESN2, SESN1, and GADD45A) and downregu- 
lation of cell cycle genes {E2F3, E2F5, DBF4, SKP2, MCM3, 
CDC27, MCM4, WEE1, MCM6, CDKN1A, PLK1, BUB1, 
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Figure 3 (a) The key connected subnetworks of overlapping signature genes of mitomycin and etoposide. Most of the genes in the maximal 
connected subnetwork are cell cycle genes, (b) Enrichment plot of the overlapping signature genes of mitomycin and etoposide. The ratio of 
enrichment is computed as -log(adjusted P value). GOBP, Gene Ontology Biological Processes annotation; KEGG, Kyoto Encyclopedia of 
Genes and Genomes. 
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Figure 4 (a) The key connected subnetworks of overlapping signature genes of monastrol and trichostatin A. (b) Enrichment plot of the 
overlapping signature genes of monastrol and trichostatin A The ratio of enrichment is computed as -log(adjusted P value). GOBP, Gene 
Ontology Biological Processes annotation; KEGG, Kyoto Encyclopedia of Genes and Genomes. 
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BUB1B, MDM2, ANAPC7, GADD45B, GADD45A, CCNA2, 
and TFDP1) in the diffuse large B-cell lymphoma LY3 cell 
lines by checking their overlapping signature genes. 

We further found that existing studies have shown that 
these two drugs have similar action mechanisms. It has 
been reported that mitomycin is activated in vivo to a Afunc- 
tional and trifunctional alkylating agent which binds to DNA 
and leads to cross-linking and inhibition of DNA synthesis 
and function. 23-25 Moreover, previous studies have shown 
that mitomycin is a cell cycle phase-nonspecific agent. Eto- 
poside inhibits DNA synthesis by forming a complex with 
topoisomerase II and DNA. This complex induces breaks 
in double-stranded DNA and prevents repair by topoisom- 
erase II binding. Accumulated breaks in DNA prevent entry 
into the mitotic phase of cell division and lead to cell death. 
The key point is that etoposide acts primarily in the G2 and 
S phases of the cell cycle. 2627 All the analysis demonstrates 
that the common action mechanism of these two drugs at 
the transcriptomic level explains why the current strategy 
can predict it well. 

The drug combination monastrol and trichostatin A, 
ranked 18 in the gold-standard list, has been predicted to 
be rank 1 9 among all (Figure 4). Through functional enrich- 
ment analysis, we found that the overlapping signature 
genes of the two drugs are involved in DNA replication and 
cell cycle pathway. Previous studies have demonstrated 
that these two drugs are the protein kinase inhibitors and 
can induce cell cycle arrest. Monastrol targets kinesin- 
related motor protein KIF11, which is required for establish- 
ing a bipolar spindle and for mitosis. Blocking of KIF1 1 can 
prevent centrosome migration and causes cell cycle arrest 
in mitosis. 28-30 Trichostatin A is a classical histone deacety- 
lase inhibitor. It has been demonstrated to induce cell cycle 
arrest, promote cell apoptosis, and inhibit metastasis. 31 We 
believe that the similar functional roles of these two drugs, 
reflected at the transcriptomic level, enable our method to 
reliably predict their synergistic effect. 

However, we found that some drug combinations cannot 
be predicted well. For example, the top drug combination on 
the gold-standard list — doxorubicin and H-7— is predicted 
to be ranked 83 by our method. In enrichment analysis, we 
found that doxorubicin mainly alters the expression of genes 
in the p53 signaling pathway, cell cycle, and DNA replica- 
tion, while H-7 mainly affects steroid biosynthesis, the spli- 
ceosome, and the ribosome. Previous studies have also 
shown that doxorubicin interacts with DNA by intercalating 
and inhibiting macromolecular biosynthesis, and induces a 
break in the DNA chain, preventing the DNA double helix 
from being resealed, thereby stopping the process of repli- 
cation, 3234 but H-7 (protein kinase C-a antibody) is a protein 
kinase inhibitor. 35 Protein kinase C has been shown to act 
as a tonic negative regulator of basal steroidogenesis in Y1 
cells by suppressing the expression of mRNA encoding the 
steroid synthetic enzymes. 36 This analysis may show that the 
predicting task based only on gene expression in treatment 
with individual drugs would fail when the effects of individual 
drugs are different. It will be necessary to incorporate more 
types of cell-response information to improve the accuracy 
of prediction. 



The drugs' different time spans to generate their maxi- 
mal effects can affect prediction 

Drugs work on cells differently, and their maximal effects can 
occur at different time points. We used a pairwise correspond- 
ing comparison between time points, which may limit the 
predicting ability. To test if the different time spans to gener- 
ate maximal drug effects affect the prediction, we added six 
cross-time point comparisons and then chose the best one for 
the ranking. Surprisingly, based on this modified strategy, we 
got a better prediction with c-index = 0.628 (Table 1), which 
is even higher than the top two methods with c-indexes of 
0.613 and 0.605 in the challenge. We found that this improve- 
ment is mainly due to the prediction of the drug combination 
blebbistatin and H-7. This combination is ranked 7 in the gold- 
standard list. In our original strategy, this combination was 
ranked 79, whereas in the modified version it was ranked 3. 
This analysis shows that considering the time point of maxi- 
mal drug effects is helpful to improve prediction accuracy. 

DISCUSSION AND CONCLUSION 

Predicting the synergistic effects of drug combinations using 
the quantitative gene expression profiles of cellular response 
to drugs is challenging. The current study shows that pre- 
diction based only on the gene expression profiles in treat- 
ment with individual drugs, not in combination, is even more 
challenging. We attempted to define the gene signatures 
that reflect the effect of all individual drugs using differential 
expression analysis and calculated the correlations of those 
gene signatures to predict the synergistic effect of drug/com- 
pound combinations. In other words, we assumed that the 
effects of individual drugs were measured in the expression 
levels of their regulated genes. The expression correlation of 
common regulated genes of two individual drugs was used to 
approximate the synergistic effect. 

The method presented here is very simple but shows very 
competitive performance compared with other methods in the 
challenge. We have shown that the differential gene signa- 
tures are more informative than randomly selected ones, and 
the method can capture common functional genes for some 
top-ranked drug pairs. The pathway analysis shows that the 
cellular response to individual drugs can partially explain how 
it works. The significant consistency between our predictions 
and the experimental gold standard demonstrates that this 
approach can effectively identify synergistic drug combina- 
tions via the cellular response to individual drugs. A potential 
limitation of this method is that it relies on the transcriptional 
signature of individual drugs and omits other complicated 
combinatorial transcriptional features of drug combinations 
in real application. The key hypothesis here is based on the 
assumption that similar effects should have similar gene sig- 
natures no matter what effects the two drugs have. However, 
this may be true for chemotherapy drugs but not true for other 
targeted-therapy drugs. Thus, the prediction strategy here 
may work for a drug combination whose individual drugs show 
related transcriptomic mechanisms but not for others. Another 
potential limitation is that the correlation-based method is not 
totally in sync with the prediction of the synergy. However, the 
performance compared to other methods demonstrates that 
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the correlation strategy can still capture some underlying con- 
nections between drugs. 

The simple method for predicting cooperative drug effects 
confirms the potential predictive ability based only on the 
quantitative cellular profiling of response to individual drugs. 
With the development of systems biology, we believe that the 
current method can be potentially improved in several aspects: 
(i) generating treated expression profiles of cell lines with drug 
combinations and building a model to compare them with those 
of each drug and (ii) increasing the number of samples and 
time points. Moreover, integrating other types of response data 
may be helpful to improve the performance. For example, Win- 
ter et al. 37 integrated transcriptional data with proteomic data 
to dissect the synergy between two multikinase inhibitors in 
chronic myelogenous leukemia cell lines. Finally, we note that 
integrating pathway and network data to improve the interpre- 
tation of synergy prediction will be necessary in future studies. 

METHODS 
Data 

We were given the genomic data on the LY3 diffuse 
large B-cell lymphoma cell line for this study. These data 
were collected both before and after treatment with iso- 
lated compounds. More detail is provided at http://www. 
the-d ream-project. org/challenges/nci-d ream-drug- 
sensitivity-prediction-challenge. 

Treated and untreated gene expression profiles were used 
in this study. The treated data included a set of 150 gene 
expression profiles corresponding to a set of 14 perturba- 
tions (14 compounds in DMSO), at three time points (6, 12, 
and 24 hours), at one concentration corresponding to the 
IC20 of the compound at 24 hours, in triplicate. The untreated 
data included the gene expression profiles of cells exposed 
to DMSO in eight replicates at all three time points. All gene 
expression profiles were profiled using an Affymetrix U219 
96 array plate (Santa Clara, CA). 

We also obtained the IC20 of the compounds and the 
synergistic activity of drug combinations assessed within 
the same experiment in triplicate to avoid differences in 
compound potency or cell line physiology in different exper- 
iments. The experimental synergistic activity was used to 
rank the computational predictions as "gold standard." 

We downloaded a human binary protein interaction net- 
work from the HPRD database (http://www.hprd.org/), which 
consists of 30,047 proteins and 41 ,327 interactions. 

Preprocessing 

We found that the gene expression profiles of the three rep- 
licates were highly correlated with each other (Pearson cor- 
relation coefficient, r > 0.99). We used the average of the 
three replicates (after normalization) as the unified profile of 
each drug. 

Method 

In the following, we describe the detailed procedure of our 
method (Figure 1b): 

1. The gene expression profile after administration of 
DMSO vehicle was regarded as the baseline profile for 



computing the differential expression after treatment 
with each drug at three different time points. 

2. For each drug, we extracted the genes with an absolute 
differential expressed z-score > 1 .96 (P < 0.05) across 
all genes for the same time point as signature genes. 

3. For each pair of drugs, we calculated their correlation 
of expression profiles based on the overlap of their sig- 
nature genes. 

4. We selected the one with maximal absolute value 
among three correlation values corresponding to three 
time points as the synergistic score for each pair of 
drug and ranked the drug pairs according to this score. 

Network-based method 

We mapped the differential expression genes onto the HPRD 
network and then extracted the maximal weighted connected 
subnetwork using the BioNet package 38 to do the same anal- 
ysis as in our correlation-based method. The detailed proce- 
dure is as follows: 

1 . Extract the significantly differentially expressed genes 
if their corresponding |z-score| > 1.96. 

2. Download a human binary protein interaction network 
from the HPRD database (http://www.hprd.org/). 

3. Assign the nodes of the significantly differentially 
expressed genes weights as their absolute z-score. 
The other nodes were assigned weights based on their 
corresponding [z-score| - 1.96. 

4. Maximal score connected subnetwork is extracted 
using the R package BioNet. 

5. Do the same downstream analysis as in our present 
method and calculate the weighted c-index. 

The c-index 

Here, a modified version of the c-index called the probabilistic 
c-index, to account for the probabilistic nature of the gold stan- 
dard, was used in authoritative scoring scripts. The probabilis- 
tic c-index was calculated by comparing the predicted rank list 
of drug pairs to the experimentally determined ranked list of 
drug combinations. Given a pair of drug combinations y'and k, 
a prediction rank was assigned a score given by the probabil- 
ity that the noisy gold standard supports this prediction. The 
probabilistic c-index is the average, for all drug pairs, of the 
probabilities that the gold standard supports the prediction. 

Functional enrichment analysis 

We adopted the DAVID online tool (http://david.abcc. ncifcri. 
gov/) to do functional enrichment analysis. The pathways and 
Gene Ontology terms with adjusted P values (Bonferroni cor- 
rection) <0.01 were selected as significant ones. 
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Study Highlights 

WHAT IS THE CURRENT KNOWLEDGE ON THE 
TOPIC? 

S Microarray data were first systematically ex- 
plored in model organisms to elucidate the drug 
mechanism of action, and coexpression analy- 
sis enables the inference of functional roles for 
genes that respond coherently to drug perturba- 
tions. Although microarray data have been ap- 
plied to explore the mechanism of drugs and the 
synergistic effects of drug combinations, there 
are currently no studies on the prediction of syn- 
ergistic effect of drug combinations using the ex- 
pression data in treatment with individual drugs. 

WHAT QUESTION DID THIS STUDY ADDRESS? 

y Can gene expression profiles in treatment with 
individual drugs be useful for predicting syner- 
gistic effects of drug combinations? 



WHAT THIS STUDY ADDS TO OUR KNOWLEDGE 

S The gene expression profiles in treatment with 
individual drugs are helpful for predicting syn- 
ergistic effects of drug combinations when the 
drug combinations show related action mecha- 
nisms at the transcriptomic level. 

HOW THIS MIGHT CHANGE CLINICAL 
PHARMACOLOGY AND THERAPEUTICS 

y This method allows elucidation of transcriptomic 
profiles for exploring synergistic effects of drug 
combinations and is applicable for determining 
target biological pathways or networks for un- 
characterized drug combinations. 
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